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Large-scale structure surveys in the coming years will measure the redshift-space power spectrum 
to unprecedented accuracy, allowing for powerful new tests of the ACDM picture as well as mea¬ 
surements of particle physics parameters such as the neutrino masses. We extend the Time-RG 
perturbative framework to redshift space, computing the power spectrum P a (k,p) in massive neu¬ 
trino cosmologies with time-dependent dark energy equations of state w(z). Time-RG is uniquely 
capable of incorporating scale-dependent growth into the P a (k, p) computation, which is important 
for massive neutrinos as well as modified gravity models. Although changes to w(z) and the neu¬ 
trino mass fraction both affect the late-time scale-dependence of the non-linear power spectrum, we 
find that the two effects depend differently on the line-of-sight angle p. Finally, we use the HACC 
N-body code to quantify errors in the perturbative calculations. For a ACDM model at redshift 
z — 1 , our procedure predicts the monopole (quadrupole) to 1% accuracy up to a wave number 
0.19/i/Mpc (0.28/i/Mpc), compared to 0.08/i/Mpc (0.07ft/Mpc) for the Kaiser approximation and 
0.19/i/Mpc (0.16h/Mpc) for the current state-of-the-art perturbation scheme. Our calculation agrees 
with the simulated redshift-space power spectrum even for neutrino masses above the current bound, 
and for rapidly-evolving dark energy equations of state, \dw/dz\ ~ 1. Along with this article, we 
make our redshift-space Time-RG implementation publicly available as the code redTime. 


I. INTRODUCTION 

Two major challenges for cosmology over the next 
decade are measuring the neutrino masses and con¬ 
straining the evolution of the dark energy density. The 
sum of the neutrino masses, a fundamental Standard 
Model parameter, is bounded from above by cosmologi¬ 
cal probes mw e m„ < 0.23 eV. Improved measure¬ 
ments of large-scale structure and the cosmic microwave 
background over the next several years will replace this 
bound with a measurement, possibly allowing us to dis¬ 
tinguish between normal and inverted neutrino mass hi¬ 
erarchies m- 

Meanwhile, searches for dark energy evolution are en¬ 
tering a decisive era. The cosmological constant A, the 
simplest model of dark energy, is completely consistent 
with current data n 0 umi]. However, a cosmologi¬ 
cal energy density which is 120 orders of magnitude be¬ 
low fundamental scales, yet coincidentally nearly equal 
to the dark matter density today, requires a great deal 
of fine-tuning iTSl - l?21 . If we consider models in which 
the dark energy density pde and its equation of state 
wde(z) = Pde/pde may vary with redshift z, the uncer¬ 
tainty in dwuE/dz is of order unity today [23]. Current 
constraints allow “early dark energy” models in which 
Pde(z) rises rapidly with z, substantially alleviating the 
tuning and coincidence problems associated with the cos¬ 
mological constant [241132] . Over the next ten to fif¬ 
teen years, constraints on divv^/dz will improve sub¬ 


stantially, allowing models with \dwEE/dz\ ~ 1 to be 
distinguished decisively from slowly-evolving equations 
of state 0 . Thus cosmology is poised to answer two fun¬ 
damental questions about the nature of the universe. 

Such powerful cosmological constraints will be made 
possible by combining probes of the expansion rate H(z ), 
including Type la supernova [33] and baryon acoustic os¬ 
cillation (BAO) 0 surveys, with measurements of the 
growth factor D(z), including redshift-space distortions 
(RSD) [331 133 and weak lensing [oH]. Here we are partic¬ 
ularly interested in the RSD, which probe the logarith¬ 
mic growth rate f(z) = —dlog(D)/dlog(l + z). These 
are measured at quasi-linear scales 10 Mpc - 100 Mpc, 
and they are minimally affected by astrophysical system- 
atics such as baryonic feedback. As a result, perturbative 
treatments of the redshift-space power spectrum are fea¬ 
sible fTITl - FTTTI . 

Time-RG perturbation theory, which directly inte¬ 
grates a system of non-linear equations for the matter 
power spectrum P{k ) (the Fourier transform of the two- 
point correlation function), was designed for models with 
scale-dependent growth factors, such as massive neutrino, 
modified gravity, and clustering dark energy models @ 33 - 
143] . It is implemented in publicly-available codes includ¬ 
ing Copter [44] and CLASS [45]. Its approach to comput¬ 
ing the power spectrum is to truncate the infinite tower 
of evolution equations for Appoint correlation functions. 
Since the continuity and Euler equations of classical fluid 
dynamics relate the time derivative of the Appoint func¬ 
tion to the (N + l)-point function, the power spectrum 
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P(k) depends upon the bispectrum B(k \, k 2 , ^ 3 ), the bis¬ 
pectrum upon the trispectrum, and so on. Time-RG 
truncates this hierarchy by neglecting the connected part 
of the trispectrum, allowing the bispectrum to describe 
the non-linear evolution of the power spectrum. 

Time-RG results were compared with those from N- 
body simulations in gU [35]. In particular, Ref. [15] 
showed that Time-RG accurately predicted the power 
spectrum even for early dark energy models which cause 
Standard Perturbation Theory to fail, and for massive 
neutrino models to which other perturbative treatments 
are inapplicable. Furthermore, Time-RG was used to test 
the approximation of Refs. S3 SB], in which neutrino 
clustering is neglected as a source for non-linear dark 
matter growth. 

In this article, we extend Time-RG to a prediction of 
the redshift-space power spectrum using the approach 
of Ref. [39], which describes higher-order corrections to 
the redshift-space power spectrum in terms of integrals 
Pb(/c,/x) and Pr(k,p,) over the bispectrum and trispec¬ 
trum, respectively. In particular, we show that PB(fc,/i) 
can be expressed as a linear combination of terms de¬ 
pending on k alone, whose time-evolution can be com¬ 
puted in the Time-RG framework. Our treatment au¬ 
tomatically includes corrections to Pb due to non-linear 
evolution, and we show that these corrections result in a 
smearing of baryon oscillations in the Pb terms as well as 
a transfer in power from larger to smaller scales. Extend¬ 
ing this calculation to massive neutrino models, we show 
that Y rn v an d the equation of state parameters affect 
the redshift-space power spectrum differently, due to the 
scale-dependent growth sourced by massive neutrinos. 

Next, we compare our calculations of the redshift- 
space power spectrum to the results of the HACC high- 
precision N-body simulations |45[ used in Ref. [15]. For 
models without massive neutrinos, we show that the 
Time-RG predictions remain accurate down to smaller 
scales than those of Ref. [35], largely due to the bet¬ 
ter small-scale behavior of Time-RG relative to closure 
perturbation theory. Finally, we demonstrate that our 
approach accurately computes the redshift-space power 
spectrum for the full range of Y rn v allowed by current 
data BED], and for models with cosmological constants 
as well as rapidly-varying dark energy. 

This article is organized as follows. Section [TT] covers 
Time-RG perturbation theory, the linear neutrino ap¬ 
proximation used here, and the basics of perturbative 
RSD calculations. Our new results are derived in Sec. ED 
and after comparing them with the calculations of [3jT we 
contrast the effects of Y m v and the equation of state pa¬ 
rameters on the redshift-space power spectrum. N-body 
simulations are used to test our results in Sec. ED Sec¬ 
tion [V] concludes that our redshift-space Time-RG calcu¬ 
lation agrees closely with N-body simulations in massive 
neutrino and evolving dark energy cosmologies, and dis¬ 
cusses possible applications. 


II. BACKGROUND 

A. Time-RG Perturbation Theory 

Consider a spatially-flat universe containing several 
non-relativistic fluids, each with density pi(r,x) and ve¬ 
locity field vi(t,x), where r and x are the conformal 
time and comoving position, respectively. caa Assume 
that the fluids interact only gravitationally. In terms of 
the density contrasts d/(r, x) = pi(r, x)/pi(r) — 1 , where 
the pi(r) are the mean densities, we can write down the 
continuity and Euler equations for each fluid as well as a 
Poisson equation coupling them: 

^ + V • V! + V • (iT/J/) = 0 (1) 

^ + Hvi + (tfj • V)Fj + V$ = 0 (2) 

V 2 <b-5fi m (T)H 2 ^/j<5j = 0 (3) 

j 

where the total cold matter density p m = Yj f jPj- Here 
is the gravitational potential and T~L = a~ 1 da/dr the 
conformal Hubble rate. In the regime of validity of cos¬ 
mological perturbation theory, the velocity fields are well- 
approximated as irrotational, V x {?,; = 0. Then each 
velocity can be described completely by its divergence 
6i = 7~L~ 1 V ■ Vi using V 2 F/ = PV0/. Henceforth we de¬ 
scribe each fluid in terms of these scalar perturbations Si 
and 9i. 

We may put the equations of motion into a more com¬ 
pact notation. Let p = log[(l + ~i n )/(l-|-z)], where Zj n 
1 is an initial redshift at which non-linear terms in the 
equations of motion are negligible. Working in redshift 
space, define <pio = exp(— 77 )dj and tpn = — exp(— p)9i. 
With primes ( , ) denoting derivatives with respect to 77 , 
and summation over repeated lower-case indices, we have 


<Pla. + ^IabPlb 

- e r, f d3 9 

J 

d3p it - 

i( 2 7r)3^c(fc, P, q) 

(4) 


xip Ib (p)ipi c (q)(2n 3 )S^(k +p + q) 

J~Uoo = —O/oi 

= 1 


(5) 

10 

= -^H m (77) 


(6) 



V W J 


H/n 

= 2 + 


(7) 

7ooi (k, q, p) 

= 7 oio (k,p, q) = {p + q)- q/i^P 2 ) 

(8) 

7111 (k, q, p) 

= {p + qfp■ 

q/(2p 2 q 2 ) 

(9) 

with all other 

'Yabc zero. 

In the first line, the 

k- 


dependence of ip on the left-hand-side has been sup¬ 
pressed, and do is the Dirac delta function. 

Next consider a universe with a single fluid. Us¬ 
ing Eq. Q, Time-RG Perturbation Theory m writes 
down the equations of motion for the power spec¬ 
trum P ab (fci)(27r) 3 (5| ) (A:i + k 2 ) = (Pa(ki)ip b (k 2 and 
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Model 

h 




n s 

^8 

Wo 

W a 

ACDM 

0.71 

0.1335 

0.02258 

0 

0.963 

0.8 

-1 

0 

yACDM 

0.71 

0.1335 

0.02258 

0.001 

0.963 

0.8 

-1 

0 

EDE 

0.6167 

0.1472 

0.02261 

0 

0.9611 

0.8778 

—0.7 

0.6722 

z/wCDM 

0.7342 

0.1543 

0.02323 

0.003 

0.8797 

0.8056 

-1.21 

-1.11 


TABLE I: List of cosmological models considered here. For 
each component J, c oj = Qj(ro)h 2 , where tq is the time today. 


bispectrum B abc (ki, k 2l /c3)(27r 3 )i5^(fci + k 2 + k 3 ) = 

(^a{ki)vb(k2)Vc{k 3 )y. 


P ab = acPcbik) - n bc P ac (k ) 

d 3 q 


Ae 71 


(27, 


hacd(.k,q,p)B bcd (k, q,p) 


+1bcd(k,q,p)B acd (k, q,p)} (10) 

P abc = ^ ad(k)B dbc Cl bd (q^ B adc Cl cd {j)^j B abd 

+2e r, [ r y a de(k, q,p)P db (q)P ec (p ) 

+7 bdeiq, P, k)P dc {p)Pea{k) 

+7 cde{p,k,q)P da (k)P eb {q)\ (11) 


where the bispectrum in the second equation is under¬ 
stood to be a function of k, q , and p: B abc = B abc (k, q,p). 

Needing to keep track of the functional dependence of 
B a bc on three different variables would make Time-RG 
perturbation theory numerically prohibitive. However, it 
is possible to recast the equations of motion in such a 
way that the time-dependent quantities depend only on 
one variable, k: 


P ab — ^ ac P bc ^ bc P ac 


47r 

(Iacd,bcd 


‘-bcd,acd 


Iacd,bef 


fq 2 dq sin ada Jacdik, q,p-)B bef {k, q,pJ) 


J (27r) 3 2 

Iacd,bef = bglacd,gef ~ &eglacd,bgf 
^/ glacd,beg d - 2e ^ A acdbe f 

1 ac d,bef _ [ d 2 dq sin ada 7 acd {k, q,p_) 

k J (2ir) 3 2 

x lbgh(k,q,P-)P ge (q)Phf(p-) 

+legh ( q,P-,k)Pgf (P-)Phb (k) 
+7 fgh(P-,k,q)P gb (k)P he {q) 


( 12 ) 

(13) 

(14) 


(15) 


Here, p± = k ± q, and the P ab , I ac d,bef , and A acd,bef 
depend on k unless otherwise specified. Since most of 
the 7 abc are zero, there are only 24 nonzero components 
of I a cd,bef a t ea ch fc, only 14 of which are independent. 

Before proceeding, we note that the treatment of non- 
linearities in Eq. (14) neglects any scale-dependence in Cl. 
In particular, in a cosmological model containing multi¬ 
ple species which cluster differently, such as cold dark 
matter and massive neutrinos, Jlio will depend on k as 


given by Eq. @. Since this is a small correction to a 
non-linear correction, it is expected to be small 155 . and 
it will be further suppressed when the density fraction of 
the second species is small. We defer discussion of the 
scale dependence of Eq. (14) and its redshift-space gen¬ 
eralization to Appendix [cr in which we show that the 
associated error is less than 1% all the way to z = 0 and 
k = 0.4 h/ Mpc even for neutrino masses m v = 0.94 eV 
much larger than current bounds. 

The Time-RG calculations of Eqs. ( I2][l5 ) may be 
sped up considerably by replacing the power spectra in¬ 
side the integral of Eq. (15) by the linear-theory power 
spectra. Since these scale in a known way with the 
growth factor D(z) and its logarithmic derivative f(z) = 
—dlog.D/dlog(l + z), the integral need not be repeated 
at each time step. This is known as the “1-loop” ap¬ 
proximation of Time-RG, since it is equivalent to stan¬ 
dard perturbation theory at the 1-loop level [H]. On a 
standard 8 -core desktop computer, 1-loop Time-RG and 
other 1-loop computations take ~ 1 minute, full Time- 
RG takes ~ 1 hour, and 2-loop perturbation theories take 
~ 1 day to compute the power spectrum over 100 k values 
in the range 0.001 h/ Mpc < k < 1 h/ Mpc. 


Figures [T] and [2] show our results for a massless- 
neutrino model with cosmological constant (ACDM) and 
an early dark energy model (EDE) whose parameters are 
given in Table |T] power spectra have been divided by the 
“no-wiggle” power spectrum P nw (k) of Ref. 52]. Time- 
RG is compared to the N-body simulations of Ref. 46], 
as well as to the 1- and 2-loop closure theory calcula¬ 
tions of Refs. [55H55] implemented in the Copter code 
[U). Full and 1-loop Time-RG are very similar in this 
range of k, with the 1-loop power spectrum slightly larger 
at k > 0.25 h/ Mpc. While 2-loop closure theory is 
highly accurate, its computation time makes it difficult 
to use in an analysis exploring a large parameter space. 
Meanwhile, 1-loop closure theory is comparable to 1- 
loop Time-RG in accuracy as well as running time in 
the k < 0.2 h/ Mpc range. Although all perturbative 
calculations break down for sufficiently large k. one ad¬ 
vantage of Time-RG is that it diverges relatively slowly 
from N-body calculations, remaining within « 10% of 
the N-body power spectrum up to k = 0.3 h/ Mpc for 
2 > 0.5 in both figures. Figure [3] a comparison to the 
high-resolution ACDM N-body simulation of Ref. 46] , 
shows that both versions of Time-RG are correct in the 
range k < 1 h/ Mpc to < 10% for z > 1 and to < 5% 
for z > 2. The figure also shows two other trends: first, 
that full and 1-loop Time-RG differ substantially only for 
k > 1 h/ Mpc; and, second, that full Time-RG underes¬ 
timates the small-scale power at high z while overesti¬ 
mating it at low z. Figure [4] compares Time-RG to the 
velocity power spectrum Pgg and the cross spectrum Pgg. 
Since the velocity field can only be determined to scales 
k ~ 0.5 h/ Mpc given our simulation resolution, we trun¬ 
cate the figure there. Evidently Time-RG is accurate at 
the 10% level to k = 0.3 h/ Mpc for both power spectra. 
It is not quite as accurate as it was for the density power 
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0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

wave number k [h/Mpc] wave number k [h/Mpc] 

FIG. 1: Power spectrum Pss for the model ACDM. (Left) Points with error bars show the N-body simulations of Ref. [323 5 an< i 
the inner (outer) shaded regions are within 5% (10%) of the simulation points. Full and 1-loop Time-RG are compared with the 
closure theory calculations of Ref. [39). Note that at the wave numbers shown, full and 1-loop Time-RG are indistinguishable 
except at z = 0 and k > 0.25 h/ Mpc. (Right) Fractional errors associated with the perturbation theories on the left. The 
shaded region shows the 2<r N-body simulation error bars. 

spectrum, but once again, Time-RG diverges from the 
N-body spectra smoothly. 


We conclude by addressing an error in the original 
Time-RG algorithm [JT]. As noted above, there are 14 in¬ 
dependent components of the A ac d t bef(k) integral at each 
k ; however, the original work included only 12 of them, 
as pointed out in Ref. [56] . The erroneous version was 
implemented in Copter, and Ref. [44] reported a discrep¬ 
ancy between simulations and the Time-RG predictions 
of the density-velocity cross-spectrum as well as the ve¬ 
locity auto-spectrum. In Figure [5] we add the corrected 
Time-RG power spectra to their figure. Including all of 
the A ac d^bef increases Pse and Pgg , improving agreement 
with simulations in the k < 0.1 h /Mpc region. This is 
significant because the velocity power spectra are impor¬ 
tant to calculations of the redshift-space power. 


B. Massive neutrinos: Linear approximation 

The coupled system of equations for multiple, non¬ 
linear, non-interacting fluids was described in Sec. m 
Since baryonic gas dynamics and feedback effects are only 
important at small scales, beyond the reach of perturba¬ 
tion theory, we regard the cold dark matter (CDM) and 
baryons as one single fluid, labelled “CB.” We also con¬ 
sider a second fluid “p” consisting of three equal-mass 
neutrino species, which behave as a warm dark matter 
component. Though a generalization to more neutrino 
species of different masses is straightforward, we focus 
on the simplest case here. Henceforth a subscript “m” 
refers to CDM, baryons, and neutrinos together, so that 
= ficDM + + fi„. The sum of neutrino masses is 

related to the neutrino density fraction oc„ = Q Vt oh 2 by 
J2 m u = 94 Wj, eV. 

Since neutrinos do not cluster strongly, they can fur¬ 
ther be approximated as a purely linear fluid in Time-RG 
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FIG. 2: Same as Figure [l] but for the early dark energy model EDE from Table [I] 


perturbation theory [41j [42|. Thus we may use a linear 
code such as CAMB or CMBFAST [777TPoTJ] to compute 
the coupling matrix Hcb- Dropping the “CB” subscript 
of we have 

fiio w -^n m (?y) [fcB + & v {r}, k)] (16) 

8&v=fv& lin,i//<Jlin,CB- (17) 

Using ^iin,CB rather than i5cb in the denominator of & v 
makes only a negligible difference [J2j. Non-linear evo¬ 
lution of is also not the dominant source of error in 
this approximation 51, 60, 61 . Thus Time-RG can in¬ 
corporate massive neutrinos, or other species leading to 
scale-dependent CB growth, without much difficulty. 

In the two-fluid case, the total matter perturbation is 
<5m = (pcb<>cb + Pv8v)/p m = /cb^cb + f v 8 v . In the lin¬ 
ear neutrino approximation, the non-linear matter power 
spectrum is given by 

P m (k) « /cb-Pcb + ZfcBfvy/PcBPlin,v + flPlin,w (18) 

Figure (el compares P C b, P\\n,v, and ^/PcB-Piin^ to the 
linear CB power spectrum. 


Based on the figure, the linear-neutrino approxima¬ 
tion should apply even beyond the quasi-linear regime 
k ~ 0.1 h/Mpc. Non-linear effects in each fluid become 
important when Aj in ~ 1, as seen by comparing the 
CB x CB and linear curves. Even for a rather large 
neutrino mass at z = 0, the neutrino auto-spectrum 
remains below 0.1, while the cross-spectrum plateaus 
around 0.25. If A 1(h) — Af inv (k) ~ Af lnv (k), then the 
leading-order fractional correction to Eq. ( fl8| ) will be 

~ {U/fcB)^A 2 n ^/A 2 CB A 2 n ^, which is at most ~ 1%. 

Moreover, non-linear corrections to A 2 should diminish 
in importance for k 3> 0.1 h/ Mpc, allowing a variant of 
this approximation to be used even in N-body simula¬ 
tions [ini m ■ 

C. Redshift-space distortions 

In a perfectly homogeneous universe, there is a pre¬ 
cise relation between the redshift 2 of an object and 
the comoving distance x{ z ) = Jq dz'/H(z') to that ob¬ 
ject. Peculiar velocities v sourced by density inhomo- 
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0 0.5 1 1.5 2 

wave number k [h/Mpc] 

FIG. 3: Same as Figure [T] but using the high-resolution 
TreePM simulation of [TB], accurate to much smaller scales. 



0 0.1 0.2 0.3 0.4 0.5 


wave number k [h/Mpc] 

FIG. 4: Density-velocity cross spectrum Pse and veloc¬ 

ity power spectrum Pee at z=l for ACDM. Inner and outer 
shaded regions correspond to 5% and 10% deviations, respec¬ 
tively, from the N-body points. 

geneities S distort this relation, perturbing x in the linc- 
of-sight direction f by an amount v • r/H (where ' de¬ 
notes a unit vector). If the homogeneous-universe rela¬ 
tion is used to identify redshift z with position x( z ) (“red- 
shift space”), then this peculiar-velocity effect shows up 
as a direction-dependent distortion of the power spec¬ 
trum. Such redshift-space distortions (RSD) enhance the 
line-of-sight power spectrum in the quasi-linear regime 
(“squashing”) and suppress it in the non-linear regime 
(“fingers of god”) [62]. 

A discussion of linear RSD as given in Ref. m is in¬ 
structive. A volume element at comoving position x and 
redshift z , at which point the velocity field takes the 



k (h/Mpc) 



FIG. 5: (a) Velocity-density cross-spectrum and (b) velocity 

auto-spectrum for the ACDM model of Ref. [44] at z = 0. 
Black points are their N-body simulation, and the solid (red) 
line is 1-loop Standard Perturbation Theory. Properly includ¬ 
ing all of the A ac d t bef integrals of Eq. ( |15| ) in the Time-RG 
calculation leads to an improved fit (solid blue line) at low k 
relative to the Time-RG results of Refs. [41] [44] (dot-dashed 
green line). Adapted from Figure 8 of Ref. 1441 . 


value v(z,x), will be assigned a redshift-space position 
x s = x + xv-x/H. For example, an object falling toward 
the observer (v ■ x < 0) will appear closer (|a? s | < |a?|) in 
redshift space than in real space. In the flat-sky approxi¬ 
mation, x ss r = constant, and the Jacobian determinant 
associated with the transformation from real- to redshift- 
space is J = \d 3 x/d 3 x s \ = (1 + d x v ■ r/'H)~ 1 . 

Since the total number of objects (or the total mass) 
in a given volume does not depend on the coordinates 
that we use, we can relate the redshift-space overdensity 
6 S to the real-space overdensity by 1 + S s = (1 + S) J. In 
Fourier space, for an irrotational velocity field v[k) oc k, 
this becomes S s = (5 + where /i£ = k ■ r. In linear 
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FIG. 6: Dimensionless power spectra A 2 (fc) = pjP(k) 

for a ui v = 0.01 (£] tn» = 0.94 eV) model at 2 = 0. The 
cross-spectrum Pcbxv ~ \JPcbPyui,v and the neutrino auto¬ 
spectrum P v « Piin.i/ use the linear-neutrino approximation. 


theory, S s 


(1 + //z|)<5, implying the power spectrum 


-P«,iin(fc,Mfe) = (1 + fp\fP\in(k). 


(19) 


This qualitatively describes line-of-sight power enhance¬ 
ment expected due to large-scale infall towards matter 
overdensities, but does not include finger-of-god effects 
at smaller scales. 

In order to extend such a treatment beyond linear 
theory, two types of corrections have been introduced: 
streaming models and higher-order corrections to the 
power spectrum [401163] . The simplest higher-order cor¬ 
rection is to replace Pi; n in Eq. (19) by the non-linear den¬ 


sity power spectrum, while the Scoccimarro anzatz |38| 
replaces all three linear power spectra with their non¬ 
linear equivalents in order to damp the redshift-space 
power spectrum: 


where Pn n is the linear matter power spectrum. The 
convention in the literature is to “absorb” the f 2 7~L 2 into 

the a 2 , so that typical values of cr.„ jlin = f Pi in dk/(67r 2 ) 
for a ACDM model at 2 ~ 1 are ~ 1 Mpc /h. 

The treatments above approximate the redshift-space 
power spectrum in terms of the real-space power spectra. 
From the form of the Jacobian it is clear that higher- 
order correlation functions ought to play a role; however, 
the gradient term in the denominator of J suggests that 
a simple Taylor expansion will be badly-behaved when 
9 ~ 1. Instead, Ref. [35) uses spatial homogeneity to 
derive an exact formula for P s without a badly-behaved 
Jacobian denominator: 


P s (k) = / d E 


ik-rAu 


6{y) + (r-V)u(y) 


5{y') + (r • V)u(y') 


(23) 


where u{y) = —v(y) ■ f/%, A u = u(y) — and 

x = y — y'- Recognizing that the Kaiser “squashing” and 
finger-of-god effects are coupled and cannot be treated 
separately, Taruya, Nishimichi, and Saito (TNS [39]) ap¬ 
ply a cumulant expansion to the expectation value in 
Eq. (23) and find a series of corrections to P s in terms 


of higher-order correlation functions. The leading-order 
corrections are obtained by neglecting all higher-order 
correlation functions except for the bispectrum and the 
disconnected parts of the trispectrum, precisely the ap¬ 
proximation used in Time-RG: 

P s (fc,/i) = F fog (fa v kfj,)[Pss(k) + 2y 2 P S g(k ) + y 4 Pgg{k) 
+P B (k,y) + P T (k,y)) (24) 


PB(k,y) 

ky 


d 3 q jiq ' 


B TNS (q 7 k - q,-k) 


(27r) 3 q . 
~B TNS (q,k,-k- q) 


(25) 


Ps{k, y%) — PgS + 2/i|-P,50 + y~Pgg. (20) 

All of these non-linear corrections essentially try to re¬ 
duce the amplitude of the linear redshift space power 
spectrum to mimic the effect of non-linear structure for¬ 
mation. 

Streaming models explicitly impose finger-of-god 
suppression by multiplying P s (k 7 y r) by a function 
F{ og (fka v yP), which typically takes a simple form: 


PT(k,y) 

k 2 y 2 


d 3 q 

Jpnf 


TTNs{q)TTTsss{k — q) (26) 


#tns(&i, k 2 , k 3 ) = Bgss(ki,k 2 , k 3 ) - y\^Bg es {ki 7 k 2 , k 3 ) 


~H z Bim(ki,k 2 ,k 3 ) 

4~ Beee(k \, k 2l k 3 ) 

TTNs(fci) = -r 3 - Pse{ki) + y~ Pgg(ki) 
k\ L 


(27) 

(28) 


Pfog(x) — 


exp(— x 2 ) 
(1 + x 2 )~ 1 


(Gaussian) 

(Lorentzian). 


( 21 ) 


Here 3a 2 , is the trace of the velocity-dispersion tensor. In 
practice, it is either fit to the data at each z, or approx¬ 
imated by its linear-theory value 


°Vlin — o 


d 3 k „ pH 2 

(2n) 3 vv 6t r 2 


dkP lin (k) (22) 


Here and henceforth, y is assumed to mean y^ unless 
otherwise labelled. A thorough study of approximations 
to P s finds this approach to be the most successful at 
matching N-body calculations and at providing an unbi¬ 
ased estimate of the growth rate / m- 

In practice, Ref. [32] computes the power spectra in 
Eq. (24) using closure theory, the bispectra in Eq. (25) 


using the tree-level approximation of Ref. [64], and the 
power spectra in Eq. (26) using linear theory. In the next 



























Section, we compute all of these terms together in Tirne- 
RG perturbation theory, making only the standard Tirne- 
RG approximation of neglecting the connected trispec¬ 
trum. There are three potential advantages to this ap¬ 
proach: 

1. Time-RG is easily-generalized to multi-component 
models with scale-dependent growth; 

2. Time-RG is well-behaved at k < 1 h /Mpc for z > 1, 
with errors < 10 %; and 

3. errors may be estimated by approximating the con¬ 
nected trispectrum, as in Ref. [56]. 

Regarding the second of these, the most promising per¬ 
turbative RSD treatments employ a streaming function 
Jfog, a phenomenological model which is not expected to 
be correct beyond the ~ 5% level at quasi-linear scales 
k > 0.1 h/ Mpc. Thus percent-level accuracy in the com¬ 
putation of the power spectra is unnecessary; much more 
important is the calculation of P(k) at the 5% — 10% 
level over a larger range of scales. 


III. REDSHIFT-SPACE DISTORTIONS IN 
TIME-RG PERTURBATION THEORY 


A. Pb and Pt in Time-RG 


Further analysis can substantially simplify the com¬ 
putation of Pb and Pp. Though we defer the details 
to Appendices [A] and |B] we summarize the results here. 
First, the /Lt-dependence of P s (fc, n) can be separated out, 
yielding a polynomial in fi 2 with /c-dependent coefficients: 

P s (k, n) = Pfog {f&vkp)[Po(k) + P 2 (k)A + P 4 (fc) A i 4 

+Pe(k)p 6 + Ps(k)A}- (29) 

Both P B and P T contribute to P 2 , P 4 , and P 6 , while Pp 
also contributes to P§: 

P B (fc,Ai) = Pb, 2{k)fi 2 + P B , 4 (fc )^ 4 + P B , 6 (fc )^ 6 (30) 

Pr(fc, l L ) = PrAk )^ 2 + P TA {k)ii 4 

+P T Ak)A + PT, S (k)A- ( 31 ) 


Second, rather than keeping track of the full func¬ 
tional dependence of the bispectrum in order to integrate 
Eq. (25) for the P B y, we can break up the bispectrum 
integrals into a series of terms Q abc , with —1 < £ < 3, 

which depend only on k and 77 . In the Time-RG frame- 

(t) 

work, the evolution equations for the Q\l c follow from 
Eq. 111]). Thus the Q^ c are analogous to the I a cd,bef{k) 
in reahspace Time-RG ( fl3| ). Appendix [A] expresses the 
P Bj as linear combinations of kQ^ c and then derives 
these evolution equations, which require the computa¬ 
tion of a two-dimensional integral analogous to A ac d,bef 



FIG. 7: Bispectrum integrals (—1 ) l kQ^ c (k) with a = b = 
c = 8, for ACDM at z = 0. Thick lines show full Time-RG, 
while thin lines show its 1-loop approximation. The other 
components are similar in magnitude as well as fc-dependence. 
The linear power spectrum is included for comparison. 


at each time 


step. Figure [?] compares kQ 


W 

abc 


to the lin¬ 


ear power spectrum, showing that non-linear contribu¬ 
tions to redslrift-space distortions are non-negligible at 
the BAO scale. The figure also shows 1-loop approxima¬ 
tions to q { 2 c . 

Meanwhile, computation of the Pt,j is straightforward. 
We describe it in Appendix [B] for completeness, but the 
only change relative to Ref. [32] is that we carry out the 
computation using the non-linear power spectra. 

The full set of evolution equations for redshift-space 
distortions in Time-RG perturbation theory is then 
Eqs. ( |12|13|14|15|A7|A8 ). Since the bispectra are small 
at early times, the Q;2 C can be initialized to zero, as 
with the I ac d,bef ■ Reference [5T] finds that the resulting 
error is approximately 1 /zi n for an initial redshift z- ln . 
and our numerical results with ^j n = 200 are consistent 
with this. Redshift-space Time-RG follows the 40 terms 
QiZ as well as 14 unique I ac d,bef , so we can expect a 
fourfold increase in computation time relative to real- 
space Time-RG. Our implementation of redshift-space 
Time-RG, redTime, is available on-line at http://www. 
hep.anl.gov/cosmology/pert.html. redTime uses the 
GNU Scientific Library [65] to evolve the equations of 
motion, and the CUBA Library ( 66 ] to compute numeri¬ 
cally the multi-dimensional integrals in Eqs. ( |15|A8|B1| ). 
The redshift-space and multipole plots in this article have 
been produced using CAMB for transfer functions and 
redTime for Time-RG calculations. 


B. Comparison to the literature 

Figures [ 8 ] and [9] compare our calculations to those of 
Ref. [32] for a ACDM cosmology. The correction terms 
Pbj and Ptj are directly compared in Fig. [ 8 j which 
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wave number k [h/Mpc] 

FIG. 8: Non-linear corrections PB,j{k) and Pr,j{k) to 

the redshift-space power spectrum. Time-RG calculations 
are shown as thick lines, and the results of [39] are thin 
lines. ACDM is shown at: (a) z = 3; (b) z = 1; (c) z = 0. 


shows that the two sets of results agree closely at early 
times and large scales. Ref. [35] predicts somewhat larger 
corrections at k < 0.2 h/Mpc, while Time-RG corrections 
become larger in magnitude around k ~ 1 h/Mpc. 


Multipole moments of the redshift-space power spec¬ 
trum P s (k,/j,) are found by projecting Eq. (24) onto a 



(b) 



FIG. 9: Ratios of the (a) quadrupole and (b) hexadecapole to 
the monopole of the redshift-space power spectrum P s (k,p), 
shown for the model ACDM. Time-RG calculations are shown 
as thick lines and TNS calculations as thin lines, with <r„ taken 
from linear theory. 


basis of Legendre polynomials for even £: 

2 /'+l f 1 

Ph,i( k ) = — 2 — J d iJ,P s {k,ij)Z?t(ji). (32) 

Defining P L ,<? = m n = f\dp p 2n F iog , 

a = fka v , and the Legendre coefficients pa such 
that S?t{x) = J2iPu x2 ‘ f° r even b we have = 

‘^ 2 ^YliPe,i rn i+j- F° r a Gaussian Ff og (ap), m ® = 
7 ( 2 ra 2 fl ,Q 2 )/a 2n+1 , where 7 ( 0 , 2 ) = e~ t t a ~ 1 dt is the 
incomplete gamma function. For a Lorentzian Pf og , 
TOq = 2 a _1 arctan(a), and for n > 1 we have the re¬ 
cursion relation a 2 mjj = 2 /( 2 n — 1 ) — m^_ 1 . 

Ratios of the quadrupole (£ = 2) and hexadecapole 
(£ = 4) to the monopole {£ = 0) are shown in Fig. [ 9 ] for 
the Time-RG and TNS [35] calculations. At early times 
and large scales, the quadrupole and hexadecapole ra¬ 
tios approach their linear Einstein-de Sitter values 50/49 
and 6/49, respectively. At smaller scales, Time-RG and 
TNS agree at early times, while TNS predicts substan¬ 
tially higher quadrupole-to-monopole and hexadecapole- 
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FIG. 10: Logarithmic contours of the redshift-space power 
spectrum P,(fcn , k±) at 2 = 1. The equation of state parame¬ 
ters wo and w a are varied from their ACDM values wo = — 1, 
w a = 0, with all other parameters held fixed. 


to-monopole ratios at late times. This is due in part to 
the fact that Time-RG overestimates the late-time non¬ 
linear power spectrum Pgg, while the closure theory cal¬ 
culation of Ref. [39| underestimates it [44] . 


C. Redshift-space power spectrum 


The two main aims of this article are to compute the 
effects of massive neutrinos and evolving dark energy 
on the redshift-space power spectrum P s (fc||, fcj_), and to 
demonstrate that our results are consistent with N-body 
simulations. We are now in a position to do the first of 
these. 

Figure [lO] shows the redshift-space power spectrum for 
a range of dark energy models. As either u> 0 or w a is in¬ 
creased, leading to an overall increase in w(z ), the power 
spectrum increases both parallel and perpendicular to the 
line of sight. Effectively the perpendicular direction con¬ 
strains the growth factor D{z) while the //-dependence of 
P s constrains the derivative / = —dlogH/dlog(l + z). 

Meanwhile, Figure [TT| varies the neutrino density frac¬ 
tion in ACDM models while keeping the early-time power 
spectrum normalization constant. The resulting varia¬ 
tion in <7g, the late-time normalization, is the largest 
contributor to the differences among the power spec¬ 
tra in Fig. 11 Since Cosmic Microwave Background 


(CMB) measurements constrain the early-time normal¬ 
ization, neutrino constraints are strongest when CMB 
and large-scale structure data are combined [5]. 

Aside from differences in cr 8 , the models in Fig. m dif¬ 
fer in their redshift-space dependence to a much greater 


FIG. 11: Logarithmic contours of the redshift-space power 
spectrum P s (k^,k±) at 2 = 1. The neutrino density pa¬ 
rameter oj v = Q, v fih 2 is varied from its ACDM value of 
zero, with the early-time (large-scale) power spectrum nor¬ 
malization equal for all ui„, and all other parameters fixed at 
their ACDM values. This normalization convention corre¬ 
sponds to as of 0.8, 0.754, 0.705, 0.660, 0.618, and 0.580 for 
of 0, 0.002, 0.004, 0.006, 0.008, and 0.01, respectively. 


extent than those in Fig. [To] Figure[l2]presents this effect 
in a different way by showing the fractional changes in 
the quadrupole caused by varying w„, wq, and w a relative 
to their massless-neutrino ACDM values. This depen¬ 
dence arises from the small-scale suppression of growth 
by massive neutrinos, which enters the linear redshift- 
space power spectrum through a scale-dependence in 
f(z,k) = —d\ogD(z,k)/d\og(l + z), a suppression of 
non-linear corrections, and a greater linear-theory veloc¬ 
ity dispersion. Thus an analysis of redshift-space distor¬ 
tions provides one more method for distinguising between 
the cosmological effects of massive neutrinos and evolving 
dark energy. 


IV. COMPARISON WITH SIMULATIONS 
A. N-body simulations with HACC 

In order to test the accuracy of our perturbative calcu¬ 
lation, we have run a suite of N-body cosmological simu¬ 
lations using the Hardware/Hybrid-Accelerated Cosmol¬ 
ogy Code (HACC), described in Ref. [35] , For each model 
in Table [T] we ran a high-resolution TreePM simulation 
with 3200 3 particles and a box size of (2.1 Gpc) 3 . Ad¬ 
ditionally, we averaged the results of 16 PM runs, each 
with 512 3 particles and a box size of (1.3 Gpc) 3 , in or¬ 
der to reduce the simulation error at quasi-linear scales 
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wave number k [h/Mpc] 



FIG. 12: Sensitivity of the quadrupole to changes in u)„, wo, 
and w a relative to their values in model ACDM at (a) z = 1 
and (b) 2 = 0. Shown for each model is the ratio of its 
quadrupole to that of ACDM, assuming linear theory for a v . 


k ~ 0.1 h/Mpc. Redsliift-space power spectra P s (k\\,k±) 
were computed in the distant-observer approximation. 
Assuming the observer to be located far from the simu¬ 
lation volume along one of the three coordinate axes, we 
used each particle’s line-of-sight velocity to shift its po¬ 
sition from real space to redshift space. Fourier transfor¬ 
mation of this shifted particle distribution yields a mea¬ 
surement of P s , which can then be averaged over multiple 
simulation runs and all three lines of sight. 

Our simulations treat massive neutrinos using the lin¬ 
ear approximation of Refs. S3 US], as described in 
Sec. |IIB| Particles in the simulation represent only the 
baryons and cold dark matter. After their non-linear 
power spectrum has been computed, the linear neu¬ 
trino power (calculated using Boltzmann integrators such 
as CAMB and CMBFAST [57H591 [67]) is added using 
Eq. (181. This approximation was shown in Ref. [35] to 
agree well with perturbative results, with a discrepancy 
at large z and u> v that can approximately be corrected by 
including the scale-dependence of the linear growth fac¬ 
tor. For an alternative approach to simulating massive 
neutrino cosmologies, see, for example, M- 



wave number k [h/Mpc] 


FIG. 13: Monopoles (thick lines) and quadrupoles (thin 
lines) of P s (k, fi) computed using Time-RG, TNS 03, non¬ 
linear Kaiser, and N-body methods, for the m v = 0 
model ACDM. 


B. m„ = 0 ACDM and a v fitting 


Perturbative calculations for the massless-neutrino 
model ACDM are compared with N-body simulations in 
Figure [13] At each redshift, the velocity dispersion pa¬ 
rameter a v in the Lorentzian streaming function is fit to 
the N-body monopole and quadrupole by minimizing % 2 . 
The N-body error in each k bin is the quadrature sum 
of the sample variance and the run-to-run standard de¬ 
viation in the 16 simulation runs, divided by \/l6. Since 
errors in the perturbative calculations in the fully non¬ 
linear regime would bias the o v fitting, resulting in a 
worse fit at low k , we restricted the fitting procedure 
to the range k < fc max = 0.2 h/Mpc. Even then, at 
2 = 0, errors in the perturbative calculations appear to 
bias a v . Time-RG in real space overestimates P{k ) at 
k ~ 0.2 h/Mpc, leading to an overestimate of a v . hence 
an underestimate in the quadrupole. Closure theory at 
1-loop, used in the TNS calculation, has the opposite 
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TABLE II: Wave number k [/i/Mpc] below which each pertur¬ 
bation theory is accurate to the given accuracy level (Acc.) 
in the ACDM massless-neutrino model. “NL Kaiser” uses 
the Time-RG real-space power spectrum in the Kaiser RSD 
formula, Eq. (|19|). 


z 

Acc. 

Time-RG 

TNS [M] 

NL Kaiser 



o 

II 

1 = 2 

O 

II 

1 = 2 

o 

II 

1 = 2 

0 

1% 

0.13 

0.09 

0.14 

0.11 

0.08 

0.12 


2% 

0.13 

0.09 

0.15 

0.11 

0.08 

0.13 


5% 

0.18 

0.11 

0.17 

0.13 

0.17 

0.13 


10% 

0.26 

0.11 

0.18 

0.14 

0.27 

0.14 

0.5 

1% 

0.15 

0.11 

0.17 

0.15 

0.11 

0.07 


2% 

0.19 

0.11 

0.17 

0.15 

0.12 

0.07 


5% 

0.21 

0.29 

0.19 

0.15 

0.32 

0.15 


10% 

0.33 

0.33 

0.23 

0.20 

1.17 

0.19 

1 

1% 

0.19 

0.28 

0.19 

0.16 

0.08 

0.07 


2% 

0.20 

0.33 

0.19 

0.20 

0.12 

0.07 


5% 

0.37 

0.40 

0.23 

0.21 

0.70 

0.19 


10% 

1.07 

0.47 

0.28 

0.32 

0.92 

0.20 

2 

1% 

0.36 

0.24 

0.21 

0.28 

0.21 

0.19 


2% 

0.43 

0.30 

0.23 

0.29 

0.26 

0.19 


5% 

0.56 

0.43 

0.29 

0.34 

0.34 

0.20 


10% 

0.72 

0.54 

0.34 

0.41 

0.49 

0.25 

3 

1% 

0.29 

0.28 

0.23 

0.28 

0.21 

0.19 


2% 

0.36 

0.31 

0.28 

0.29 

0.26 

0.19 


5% 

0.50 

0.40 

0.34 

0.36 

0.34 

0.23 


10% 

0.66 

0.51 

0.41 

0.43 

0.49 

0.29 


effect. For z > 1, both perturbative calculations agree 
well with the simulation. The cr„ fitting procedure is 
discussed further in Appendix 0 which shows that our 
results are insensitive to the value of & max within the 
range 0.15 /i/Mpc < fc max < 0.25 h/ Mpc. 

Tabic |TT| uses the data from Fig.[l3]to list the maximum 
k up to which each perturbative calculation agrees with 
simulations to a given accuracy level. The results follow 
the broad patterns expected from Fig. [l] for percent- 
level accuracy at 0 < 0.5, the closure-theory-based ap¬ 
proach of TNS is somewhat better than Time-RG, while 
for 5%-10% accuracy levels and at higher z, Time-RG is 
preferred. In particular, we note that Time-RG diverges 
from the N-body power spectra rather smoothly, with 
no model-dependent catastrophic errors which could bias 
data analyses or forecasts. Meanwhile, the non-linear 
Kaiser model, which uses the Time-RG non-linear power 
spectrum instead of the linear one in Eq. , fails to 
match accurately the N-body quadrupole even at z ~ 1. 

Similar patterns can be seen in Figure [14} which shows 
the two-dimensional power spectrum P s (fc| , k±) for simu¬ 
lations, Time-RG, TNS, and the non-linear Kaiser model. 
The trends seen in Fig. [13] are evident here: Time-RG 
overestimates the monopole, particularly at z = 0; TNS 
underestimates the monopole; and the Kaiser model sub- 



0.05 0.1 0.15 0.2 0.25 0.3 


(b) 


o 

CL 





0.05 0.1 0.15 0.2 0.25 0.3 

k|| [h/Mpc] 


FIG. 14: Logarithmic contours of the redshift-space power 
spectrum P s (fcn, k±) for the massless-neutrino model ACDM, 
at (a) z = 0; (b) z = 1. Points show N-body calculations, 
thick dashed (colored) lines show Time-RG, solid black lines 
show TNS, and thin dashed black lines show the non-linear 
Kaiser approximation. 


stantially underestimates the quadrupole, even at z = 1. 


C. Evolving dark energy and massive neutrinos 

Finally, we apply our calculation to the last three mod¬ 
els in Table |Tj an early dark energy, a massive-neutrino 
ACDM model, and a dark energy with w a = —1.11 and 
massive neutrinos. Since next-generation galaxy surveys 
will go beyond the BOSS redshift [2] of z = 0.57 to z ~ 1, 
we focus here on z = 1 and 0. Evidently from Figure [15] 
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FIG. 15: Monopoles (thick lines) and quadrupoles (thin lines) of the redshift-space power spectra at (left) z = 1 and 
(right) z = 0 for (a) EDE, early dark energy with massless neutrinos; (b) i^ACDM, with Y2, m v = 0.094 eV; (c) vwCDM, a 
rapidly-evolving dark energy with Yh = 0.29 eV. In each case, points represent N-body calculations, solid (red) lines show 
the Time-RG calculation, and dashed (green) lines show the non-linear Kaiser model. 
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TABLE III: Wave number k [/i/Mpc] below which each perturbation theory is accurate to the given accuracy level (Acc.) for 
the models in Table [I “NL Kaiser” uses the Time-RG real-space power spectrum in the Kaiser RSD formula, Eq. (|19|) . 


EDE 


i/AGDM 


i/wCDM 


z 

Acc. 

Time-RG 

NL Kaiser 

Time-RG 

NL Kaiser 

Time-RG 

NL Kaiser 



II 

O 

1 = 2 

O 

II 

1 = 2 

O 

II 

<M 

II 

O 

II 

1 = 2 

O 

II 

1 = 2 

O 

II 

1 = 2 

0 

1% 

0.11 

0.06 

0.04 

0.10 

0.11 

0.06 

0.05 

0.12 

0.12 

0.06 

0.06 

0.11 


2% 

0.12 

0.07 

0.04 

0.10 

0.13 

0.06 

0.06 

0.12 

0.13 

0.06 

0.06 

0.11 


5% 

0.16 

0.07 

0.12 

0.10 

0.16 

0.06 

0.16 

0.13 

0.15 

0.08 

0.15 

0.11 


10% 

0.20 

0.07 

0.21 

0.11 

0.25 

0.11 

0.27 

0.14 

0.25 

0.11 

0.26 

0.13 

0.5 

1% 

0.16 

0.10 

0.07 

0.15 

0.13 

0.06 

0.05 

0.06 

0.14 

0.08 

0.04 

0.04 


2% 

0.17 

0.10 

0.07 

0.15 

0.15 

0.11 

0.05 

0.06 

0.15 

0.08 

0.04 

0.06 


5% 

0.22 

0.10 

0.22 

0.16 

0.21 

0.12 

0.32 

0.14 

0.21 

0.11 

0.35 

0.06 


10% 

0.25 

0.15 

0.33 

0.16 

0.34 

0.49 

1.17 

0.17 

0.33 

0.39 

1.12 

0.16 

1 

1% 

0.18 

0.10 

0.07 

0.09 

0.19 

0.17 

0.06 

0.06 

0.15 

0.13 

0.06 

0.06 


2% 

0.18 

0.10 

0.07 

0.09 

0.19 

0.17 

0.06 

0.06 

0.21 

0.17 

0.06 

0.06 


5% 

0.25 

0.16 

0.47 

0.18 

0.36 

0.45 

0.70 

0.15 

0.89 

0.46 

0.52 

0.16 


10% 

0.38 

0.39 

1.20 

0.19 

1.07 

0.52 

0.92 

0.20 

1.00 

0.52 

0.75 

0.19 

2 

1% 

0.25 

0.41 

0.12 

0.09 

0.45 

0.30 

0.06 

0.08 

0.37 

0.25 

0.12 

0.08 


2% 

0.25 

0.55 

0.36 

0.18 

0.49 

0.35 

0.26 

0.08 

0.44 

0.31 

0.27 

0.13 


5% 

1.06 

0.83 

0.53 

0.19 

0.62 

0.47 

0.36 

0.20 

0.54 

0.37 

0.35 

0.21 


10% 

1.13 

1.10 

0.85 

0.25 

0.79 

0.62 

0.53 

0.26 

0.69 

0.54 

0.49 

0.26 

3 

1% 

0.62 

0.39 

0.26 

0.19 

0.36 

0.34 

0.26 

0.20 

0.37 

0.37 

0.23 

0.13 


2 % 

0.69 

0.44 

0.33 

0.19 

0.45 

0.38 

0.26 

0.20 

0.44 

0.37 

0.28 

0.21 


5% 

0.86 

0.58 

0.44 

0.25 

0.60 

0.51 

0.38 

0.25 

0.58 

0.52 

0.40 

0.26 


10% 

1.05 

0.75 

0.65 

0.31 

0.77 

0.65 

0.55 

0.32 

0.74 

0.65 

0.57 

0.33 


the Time-RG RSD calculation presented here agrees well 
with N-body calculations of the monopole, quadrupole, 
and two-dimensional redshift-space power spectrum for 
a wide range of models at z = 1. By contrast, the Kaiser 
RSD model, Eq. (191 applied to the non-linear real- 
space power spectrum, substantially underpredicts the 
quadrupole at quasi-linear scales. At 2 = 0, neither per¬ 
turbative calculation is accurate beyond k ~ 0.1 /i/Mpc 
for any of the models. Table III lists the accuracy of 
both perturbative methods for all three models over the 
redshift range 0 < z < 3. 


V. CONCLUSIONS 

Cosmological surveys over the next decade will mea¬ 
sure the sum of neutrino masses, and should either de¬ 
tect or decisively exclude order-unity variations in the 
dark energy equation of state. Observations are made in 
redshift space, with the redshift dependent on the line- 
of-sight velocity of an object as well as its distance. Since 
peculiar velocities are sourced by overdensities, redshift- 
space distortions provide additional information about 
the scale-dependent growth of large-scale structure. Per¬ 
turbative techniques for calculating the redshift-space 
power spectrum in models with scale-independent growth 
were reviewed in Ref. }40j , with the method by Taruya, 


Nishimichi, and Saito [US] (TNS) proving the most effec¬ 
tive. 


In this article we have extended the TNS approach 
to models with scale-dependent growth using the Time- 
RG perturbation theory, in which higher-order contribu¬ 
tions to the power spectrum are described in terms of 
the bispectrum. We have decomposed the higher-order 


corrections found by TNS into integrals Q^ c (k) over the 
bispectrum. Using the Time-RG framework, we have de- 
rived the evolution equations for the Q^ bc and computed 
the redshift-space power spectrum in models with mas¬ 
sive neutrinos as well as rapidly-evolving dark energy. 
Finally, we have confirmed the accuracy of our calcu¬ 
lations by comparing them to N-body simulations con¬ 
ducted using the HACC code. We compare our results 
to TNS and simulations in Figures [l3| and |Td| as well as 
in Table [TT} Figure [l5| and Table III show that our results 
agree closely with N-body simulations for a wide range of 
models with massive neutrinos and rapidly-evolving dark 
energy equations of state. 


Our work is applicable in a variety of ways. Since our 
results are accurate at the 10% level over a fairly large 
range of wave numbers, particularly for z > 1, they can 
be used to forecast constraints from large-scale structure 
surveys. Forecasts based upon linear perturbation theory 
are truncated at k ~ 0.1 h/ Mpc, since linear theory at 
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smaller scales overestimates the amount of information 
available from the BAO peak [3]- Accuracy can be im¬ 
proved by building an emulator combining our perturba¬ 
tive calculation at large scales with interpolated N-body 
power spectra at small scales (5SJ1M]- Finally, perturba¬ 
tive techniques allowing for scale-dependent growth are 
applicable to modified gravity as well as to massive neu¬ 
trinos. 


Appendix A: Evolution of Pb in Time-RG 


Rather than computing the full functional depen¬ 
dence of the bispectrum B(k,q,p) on three different 
wave numbers, we may decompose PpJk, //) into a set of 
fc-dependent bispectrum integrals whose time-evolution 
may be computed within the Time-RG framework. Ex¬ 


panding the expression (251, we have 


Pb 

kp. 


d 3 q p$ 
(27r) 3 q L 


Bses(k, q,p+) 


—p 2 Bggg(k, q,p+) - pl + B 5 gg(k,q,p+) 
+p 2 p,p + Bggg(k,q,p+) - \p+ p_] (Al) 


(£) 

where the Pbj and Q abc are understood to depend on k. 

Next, we study the evolution of the Q abc in Time-RG. 
Since the evolution of the B abc (k,q,p) only depends on 
the other components of B abc (k,q,p), the evolution of 
the Q abc for fixed i will only depend on the other com- 
ponents of Q y abc - We can easily multiply the bispectrum 
evolution equation ( |TT] ) by (k/p dz )(q/k) cre ^\t\(cosa), in¬ 
tegrate over d 3 q , and pull the time derivative outside of 
the integral. Therefore, 


dgQal =-n ad Q ( d l-ttbdQai- n cd Q w 


I o p v dW 

abd ^ 11 abc 


(A7) 


= 


'q 2 dq sin ada k 


P\ 


(|) ^iq(cosa) 


(2t r) 3 

7 ade(k,q,p + )P db (q)P ec (p + ) 

+7 bde {q,p+,k)P dc (p + )P ea {k) (A8) 

+7 cde(p+, k,q)P da {k)P eb {q) + (-1 Y[p+ O jp_] 


where the Q^ c and R^. 


at each rj depend only on k. 


Appendix B: Computation of Pr 


where p± = k ± q implies p±pi ± = k 2 p 2 + q 2 pi±2kqpp^. 
The three-dimensional integral over q can be evaluated 
as the integral over the magnitude q , the angle a between 
q and k, and the angle /3 of q — (q ■ k)k in the plane per¬ 
pendicular to k. Since B abc (ki , k- 2 , ^ 3 ) is invariant under 
rotation in ki~ space, the only /3-dependent quantities in 
the integral are factors of raised to integer powers; 
thus, the integral over /3 can be evaluated trivially. 

Define the quantities 


Beginning with Eq. (26), we find Pt by integrating 


over the angle /3 of q in the plane perpendicular to k. 
Once again defining a as the angle between k and q, and 
defining r ab = P a fi(q)P bt g(p _), 

Ptj f q 2 dqs a da ^ . 

-^ 2 = J (2 ^) 3 p 2 Z. S -^ab for 3 = 2,4, 6,8 (Bl) 


where the S ab j are given by: 


Q l lW = 


q 2 dqsm(a)da k f q\ ae 


(2 7 T ) 3 


p\ 


(|) ^iq(cosa) 


a e = 


xB abc (k,q 7 p + ) + {-\Y[p+ Rp_] (A2) 
f sign(£) if t is odd 
[ 0 if l is even 


(A3) 


for l = — 1 , 0, 1, 2, and 3, where the are Legendre 
polynomials. Then 


d^B,2 _ oa>(~ 1 ) o/'F 1 ) ^< 7 ( 0 ) 

~ A ^S9S ^595 3 ^696 ^596 

4--O^ 2 ) — - 0 (°) + — -O^ (A41 

T.^S99 3 ^<5 99 ^^599 ^599 

Bp,A _ o/+ — L o/n(4) ^/n( 2 ) 

Jr j t — ^995 ^998 ^99S ^995 

Qsgg "* — 4(2^00 — 2Q^00 + -Q000 

4 . ®n( 3 ) _ -n (1) (as) 

g^ggg a- ^Wggg ^Wggg 
= ^ 2 Q { 9e9 ) ~ 4 Q 000 — 2 Qggg (A6) 


Sss, 2 = -4 (B2) 

Sse ,2 = -3s* 4/(44.) (B3) 

Sgs, 2 = -3s®/4 (B4) 

See, 2 = -5s® /8 (B5) 

Sss ,4 = 2 C a k/q + s 2 - 2c 2 (B6) 

Ssg ,4 = — 3s 2 + 9c Q s 2 + s 2 (s 2 — 4c 2 )|^r (B7) 
Sgs ,4 = 3c a s 2 k/q + 3s 2 (s 2 - 4c 2 )/2 (B8) 

= -4g + 44^ + 4(4 - 6c 2 )|^ (b 9 ) 

Ss 9,6 = C 3 a /qp 2 _ + 3(s 2 - 2 c 2 a )kq/p 2 _ 

+6(44- 4/3 -4/8)64/p 2 - (BIO) 

Sgs,6 = (2c 2 - 3s 2 )k/q + 6(c 2 s 2 - c^/3 - s^/8) (Bll) 


See, 6 = 3 c a s 2 a k 3 /(qp 2 _) + 9s 2 (s 2 - 4c 2 )k 2 /(4p 2 _) 
+30(c 3 s 2 - c a s*/4)kq/p 2 _ 

+45s 2 (c 2 s 2 - 24/3 - 4/12)4/(24.) (B12) 
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wave number k [h/Mpc] 


FIG. 16: Error associated with neglecting the fc-dependence 
of fli0; Eq. |6|, in the non-linear Time-RG corrections 
Eqs. ( |14|A7| >. Differences between power spectra computed 
with the maximum and minimum values of fiio are shown 
at z = 0 for a model with uj v = 0.01, corresponding to 
Tn„ = 0.94 eV. 


See,8 


(24-3^4)^ + 18(44-K- 
+ 6(4 — 544 + 

-(24 -1544+ f 44 - 14 )^ 


8 b a> v t 


(B13) 


We have used the shorthand notation s a = sin(a) and 
c a = cos(a). All components of S a b,j not listed above 
are zero. 


Appendix C: Error bounds and a v fits 
1. Scale-dependence of ft 


In the general multi-species case as well as in the lin¬ 


ear neutrino approximation of Sec. IIB the linear evolu¬ 
tion matrix ft is scale-dependent. In particular, Eq. © 
shows that the fractional change in fiio is of order f v 
in massive neutrino models. This scale-dependence has 
been neglected in real-space and redshift-space Time-RG, 
Eqs. ( |14|A7D , in which ft has been pulled outside the in¬ 
tegrals over wave number. Here we estimate the error 
associated with this approximation and show that it is 
small enough to be negligible even for neutrino masses 
several times current bounds. 

In order to place an upper bound on this fc-dependent 
ft error, we study a model similar to ^ACDM but with 
ten times the neutrino content, = 0 . 01 , corresponding 
to = 0.94 eV. The function l~2 10 (;r, k ) drops from a 

maximum value at large scales, where neutrinos cluster 
like cold matter, to a small-scale value that is smaller by 
a fraction ss /„/( 1 — fv)- Moreover, we expect the error 
to be smaller at higher z, since it only affects non-linear 
correction terms, which grow with time. 



redshift z 


FIG. 17: Best-fit velocity dispersion a v for the models in Ta- 
blc[T| fitting to a suite of 15 Particle-Mesh N-body simulations 
at scales k < fc m ax = 0.2 h/ Mpc. 


We estimate the ^-dependent ft error by computing 
power spectra using either the maximum or the mini¬ 
mum values of fl 10 in Eqs. ( 14|A7 ). Our estimate, half 
the difference between the two, is shown in Fig. [T 6 |fo r the 
power spectrum components Pj ( k ) defined in Eq. 29 For 
the dominant contributors to the redshift-space power 
spectrum, P 0 , P 2 , and P 4 , the error is less than 0.5% for 
k < 0.15 h /Mpc and less than 0.9% for k < 0.4 h/ Mpc. 
Errors in the smaller terms Pq and Ps are less than 1% 
for k < 0.25 h/ Mpc. At higher z, and for more realistic 
/„, we expect these errors to be several times smaller, 
meaning that they are negligible compared with the er¬ 
rors listed in Tables El and IIII1 


2. Best-fit a v 


Our results in Tables [XT] and III are based on fitting 
the velocity dispersion a v in the Lorentzian streaming 
function F{ 0 g(fa v k/j,) to the N-body simulations over the 
range k < fc max = 0.2 h/ Mpc. The resulting velocity 
dispersions are shown in Fig. [lT] for the models in Ta¬ 
ble |T[ Velocity dispersions for models ACDM, zvACDM, 
and zvwCDM are similar, while the early dark energy 
model EDE, which enhances structure growth, has a 
larger velocity dispersion. Also shown, for ACDM, are 
the velocity dispersion in the TNS calculation and that of 
linear theory, Eq. (22). For ACDM, both non-linear cal¬ 
culations find a v below linear theory at high z and above 
linear theory at low z. The best-fit a v values appear to 
be biased by perturbation theory errors at low z. Time- 
RG, which overpredicts late-time small-scale power, has 
a higher a v than TNS, which underpredicts it. 
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TABLE IV: Sensitivity to fc max . Wave numbers k [/i/Mpc] 
below which Time-RG perturbation theory is accurate to the 
given accuracy level (Acc.) for model I'ACDM are shown. In 
each case, the velocity dispersion <j v in the streaming func¬ 
tion F{ og (fa v kfj.) was fit to the N-body power spectrum over 
the range k < k„ mx . For comparison, Tables [TT] and [Till used 

fcmax = 0.2 h/ Mpc. 


z 

Acc. 

^max 

t —1 <M 

O II 

II ^ 

^max 

= 0.15 

1 = 2 

^max 

= 0.25 

1 = 2 

/i/Mpc 

£ = 0 

/i/Mpc 

£ = o 

h/ Mpc 

£ = o 

0 

1 % 

0.09 

0.06 

0.11 

0.06 

0.13 

0.06 


2 % 

0.11 

0.06 

0.13 

0.06 

0.13 

0.06 


5% 

0.14 

0.11 

0.15 

0.11 

0.19 

0.06 


10 % 

0.21 

0.17 

0.21 

0.11 

0.27 

0.11 

0.5 

1 % 

0.13 

0.11 

0.13 

0.11 

0.14 

0.06 


2 % 

0.13 

0.12 

0.13 

0.11 

0.19 

0.06 


5% 

0.19 

0.17 

0.21 

0.17 

0.26 

0.11 


10 % 

0.30 

0.36 

0.32 

0.40 

0.42 

0.17 

1 

1 % 

0.16 

0.17 

0.17 

0.17 

0.19 

0.14 


2 % 

0.19 

0.17 

0.19 

0.17 

0.21 

0.17 


5% 

0.32 

0.36 

0.32 

0.44 

0.96 

0.17 


10 % 

1.09 

0.45 

1.07 

0.49 

1.04 

0.59 

2 

1 % 

0.19 

0.19 

0.41 

0.17 

0.41 

0.29 


2 % 

0.28 

0.21 

0.49 

0.30 

0.49 

0.30 


5% 

0.78 

0.86 

0.60 

0.38 

0.61 

0.38 


10 % 

0.96 

0.99 

0.77 

0.58 

0.77 

0.59 

3 

1 % 

0.19 

0.18 

0.36 

0.34 

0.35 

0.34 


2 % 

0.28 

0.19 

0.45 

0.38 

0.45 

0.35 


5% 

1.08 

0.28 

0.58 

0.49 

0.58 

0.49 


10 % 

1.14 

0.53 

0.76 

0.64 

0.76 

0.62 


3. Sensitivity to fc max 


In the quasi-linear regime, errors in perturbation the¬ 
ory increase with increasing k, while uncertainties in 


N-body simulations and actual data decrease due to a 
greater number of modes. Thus one might worry that our 
a v fitting procedure is dominated by wave numbers near 
tax, which may be at the edge of the regime of validity 
of perturbation theory. Table IV tests the sensitivity of 
our redshift-space Time-RG perturbative calculation to 
fcmax for model zzACDM. The table shows wave numbers 
up to which perturbation theory is accurate for several 
accuracy thresholds, for /c max = 0.1 h/ M pc, 0 .15 h/ Mpc, 
and 0.25 h/ Mpc. Comparison to Table |III| shows that 
fcmax = 0.15 ft./Mpc, 0.2 /i/Mpc, and 0.25 ft./Mpc are 
similar. At the 1 — 2% accuracy levels and at z > 2, 
fc max = 0.1 h/ Mpc leads to a significantly worse fit. 
This is not surprising, since the N-body error bars in 
Fig. 15 (b) are large for k < 0.1 h/ Mpc, and since not 
much non-linear information is available at such large 
scales. We conclude that our results are stable over the 
range 0.15 h/Mpc < k m ax < 0.25 h/Mpc. 
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